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Chapter  1 

Flexible  Multibody  Systems 


The  Air  Force  remains  active  in  space  systems,  and  hardware,  such  as  space 
robots,  rotorcraft  and  spacecraft.  They  consist  of  subsystems  which  can  be  de¬ 
scribed  as  flexible  multibodies.  The  dynamics  and  control  of  flexible  multibody 
systems  has  been  of  interest  for  many  years  (Refs.  1-16).  Identifying,  modeling 
and  controlling  such  systems  using  various  theories  with  confidence  has  also  be¬ 
come  an  important  issue.  At  present,  a  real  need  exists  for  the  validation  and 
comparison  of  various  modeling  and  control  theories  based  on  an  actual  hard¬ 
ware  experiment.  However,  compared  to  theoretical  developments  and  number 
of  computer  programs  available,  experimental  verification  has  never  been  con¬ 
ducted.  To  this  end,  the  Phillips  Laboratory  at  Kirtland  AFB  has  constructed 
a  flexible  multibody  structure  called  PACE(Planar  Articulating  Controls  Ex¬ 
periment)  which  consists  of  2  flexible  beams  connected  in  series  with  motors  at 
both  the  hub  and  the  elbow  joint  (Fig.  1.1). 

In  recent  years,  the  interest  has  broadened  so  as  to  include  maneuvering  and 
control.  This  narrows  the  choice  of  formulations  significantly,  as  the  formulation 
must  be  consistent  with  the  control  task.  A  set  of  equations  of  motion  suited 
for  the  control  task  can  be  formulated  by  means  of  Lagrange’s  equations  for 
flexible  bodies  in  terms  of  quasi-coordinates.  The  advantage  of  this  approach  is 
that  it  yields  equations  in  terms  of  body  axes,  which  are  the  same  axes  as  those 
used  to  express  control  forces  and  torques.  In  using  the  approach  of  Ref.  12 
to  derive  equations  of  motion  for  a  chain  of  flexible  multibody  systems,  it  is 
convenient  to  adopt  a  kinematical  procedure  permitting  the  expression  of  the 
velocity  vector  of  a  nominal  point  in  a  typical  body  in  terms  of  the  velocity 
vector  of  the  preceding  body  in  the  chain.  The  resulting  differential  equations 
are  nonlinear  and  hybrid  (Ref.  13),  where  the  term  ‘‘hybrid”  implies  that  the 
equations  for  the  rigid-body  translations  and  rotations  are  ordinary  differential 
equations  and  those  for  the  elastic  motions  are  partial  differential  equations. 
Because  maneuvering  and  control  design  in  terms  of  hybrid  equations  is  not 
feasible,  the  partial  differential  equations  must  be  transformed  into  sets  of  ordi- 
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Figure  1.1:  PACE  Test  Article 


nary  differential  equations  by  means  of  a  discretization-in-space  procedure,  such 
as  the  finite  element  method  (Ref,  14)  or  a  Rayleigh-Ritz  based  substructure 
synthesis  (Ref.  15).  The  resulting  formulation  consists  of  a  high-order  set  of 
nonlinear  ordinary  differential  equations.  A  common  approach  to  control  de¬ 
sign  requires  the  solution  of  a  two-point  boundary-value  problem,  which  is  not 
feasible  for  high-order  systems,  so  that  a  different  approach  is  advisable. 

The  nonlinearity  enters  into  the  differential  equations  through  the  rigid- 
body  motions.  Indeed,  the  elastic  motions  tend  to  be  small.  In  view  of  this,  it 
appears  natural  to  conceive  of  a  perturbation  approach  whereby  the  rigid-body 
motions  can  be  regarded  as  being  of  zero-order  in  magnitude  and  the  elastic 
motions  as  being  of  first-order  in  magnitude.  This  approach  permits  dividing 
the  problem  into  a  low-dimensional  set  of  nonlinear  zero-order  equations  for  the 
rigid-body  motions  and  a  high-dimensional  set  of  linear  first-order  equations  for 
the  elastic  motions  and  the  perturbations  in  the  rigid-body  motion,  where  the 
order  is  to  be  taken  in  a  perturbation  sense.  Note  that,  because  the  zero-order 
solution  enters  into  the  first-order  equations  as  a  known  function  of  time,  the 
first-order  equations  represent  a  time-varying  system.  Moreover,  the  system  is 
subjected  to  persistent  disturbances.  The  perturbation  approach  just  described, 
first  proposed  in  Ref.  16,  was  used  in  Refs.  17-21  to  maneuver  and  control  flexible 
spacecraft. 

The  kinematical  synthesis  of  Refs.  16-21  works  quite  well  in  the  case  in  which 
the  number  of  bodies  in  the  chain  is  relatively  small.  When  the  number  of  bodies 
is  larger  than  three,  difficulties  can  be  expected,  so  that  a  different  approach  is 
advisable.  In  Refs.  14-19,  the  kinematical  synthesis  was  implemented  before  the 
derivation  of  the  equations  of  motion  was  carried  out.  In  this  paper,  we  consider 
a  procedure  whereby  the  equations  of  motion  are  derived  first  for  each  individual 
flexible  body.  Then,  the  sets  of  equations  for  the  individual  bodies  are  assembled 
into  a  global  set  by  invoking  the  kinematical  relations  described  above.  In  the 
process,  the  redundant  coordinates  and  velocities  resulting  from  considering 
the  individual  bodies  separately  are  eliminated.  It  is  convenient  to  carry  out 
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the  kinematical  synthesis  on  the  zero-order  problem  and  first-order  problem 
separately.  Implementation  of  the  kinematical  synthesis  is  based  on  recursive 
relations  that  lend  themselves  to  ready  computer  coding.  The  resulting  zero- 
and  first-order  global  sets  of  equations  are  particularly  suited  for  maneuvering 
and  control  design,  respectively. 

The  zero-order  nonlinear  equations  govern  the  maneuver  as  if  the  system 
consisted  of  articulated  rigid  bodies  where  the  maneuver  amounts  to  driving 
the  system  from  an  initial  state  to  a  final  state.  The  equations  can  be  solved 
open-loop  or  closed-loop.  For  minimum-time  maneuvering,  the  control  law  is 
bang-bang.  The  simplest  approach  is  to  carry  out  the  maneuver  by  means  of 
actuators  that  impart  predetermined  motions  to  the  substructures  relative  to 
one  another.  The  first-order  equations  govern  the  elastic  vibrations  and  the 
perturbations  in  the  rigid-body  motions.  They  contain  the  zero-order  solution 
as  a  known  function  of  time.  As  a  result,  the  system  is  time-varying.  Moreover, 
it  is  subjected  to  persistent  disturbances  caused  by  the  maneuver.  The  process 
can  be  likened  to  that  in  which  the  system  must  follow  a  reference  state.  In 
this  case,  the  reference  state  is  defined  by  the  rigid-body  maneuvering,  which 
is  characterized  by  zero  elastic  states.  Then,  the  first-order  equations  are  sim¬ 
ply  the  equations  in  terms  of  the  difference  between  the  actual  states  and  the 
reference  states,  where  this  difference  can  be  identified  as  perturbations  in  the 
state  variables.  Control  of  the  first-order  system  is  carried  out  closed-loop  and 
includes  a  disturbance  accommodation  procedure.  The  latter  depends  on  the 
type  of  disturbance,  which  in  turn  depends  on  the  type  of  maneuver  performed. 

As  a  numerical  example,  PACE  with  direct  DC  motors  is  considered,  in 
which  bang-bang  control  voltages  are  applied  to  shoulder  and  elbow  motors. 
Numerical  results  show  the  effectiveness  of  the  new  algorithm  developed  in  this 
paper. 
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Chapter  2 

Modeling  of  Flexible 
Multibody  Systems 


2.1  Equation  of  Motion 


Our  interest  lies  in  deriving  equations  of  motion  capable  of  describing  the  prob¬ 
lem  of  maneuvering  and  control  of  structures  in  the  form  of  an  articulated 
chain  of  2  flexible  substructures.  To  this  end,  we  propose  to  use  a  kinematical 
synthesis  procedure  whereby  the  motion  is  referred  to  sets  of  local  body  axes 
embedded  in  the  undeformed  substructures  and  the  motion  of  one  substructure 
is  described  in  terms  of  the  motion  of  the  preceding  substructure  in  the  chain. 
This  kinematical  procedure  obviates  the  need  for  constraint  equations. 

The  approach  used  in  this  report  is  to  derive  equations  of  motion  for  the 
individual  substructures  separately  and  then  impose  kinematical  relations  of 
the  type  described  earlier  to  obtain  system’s  equations  of  motion.  Although  the 
approach  is  used  for  the  case  of  a  two-link  flexible  body  system,  the  approach 
can  be  extend  to  the  case  of  arbitrary  N  flexible  multibody  systems. 

Let  us  consider  a  typical  flexible  substructure  moving  on  a  horizontal  sur¬ 
face  (Fig.  2.1)  and  introduce  the  inertial  axes  XY  with  the  origin  at  O  and  a 
set  of  body  axes  with  the  origin  at  S  and  embedded  in  the  undeformed 
substructure.  Then,  we  can  write  the  position  vector  of  a  typical  point  in  the 
substructure  with  the  spatial  coordinates  given  symbolically  as  follows: 

=  n’’  (ii,  +  Cj  u,)  =  h'^W,  (2.1) 

where  h  =  [I  represents  the  column  matrix  consisting  of  the  unit  vectors 
in  X-  and  Y-directions  corresponding  to  inertial  coordinates,  C,  =  0(0,)  is  the 
matrix  of  the  direction  cosines  which  is  given  by 


cOf  s9, 

—s6}  cdf 


(2.2) 
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X 

Figure  2.1:  Flexible  Body 

In  addition,  Rs  is  the  radius  vector  from  /  to  S,  and  Ug  includes  the  radius 
vector  from  5  to  a  typical  point  in  s  and  the  elastic  displacement  vector  of  the 
same  point  relative  to  the  body  axes  x^Vsy  respectively.  Thus,  we  can  write 

where  Vs  represents  the  elastic  displacement.  Note  that  Rs  is  in  terms  of  inertial 
coordinates  and  Xs  and  Vs  a,re  in  terms  of  components  along  the  body  axes.  We 
propose  to  derive  the  equations  of  motion  by  the  approach  of  Ref.  1.  In  fact, 
the  equations  of  motion  for  the  individual  bodies  are  identical  to  the  hybrid 
equations  derived  here,  so  that  here  we  merely  present  the  pertinent  material. 

The  Lagrangian  formulation  requires  the  kinetic  energy,  which  in  turn  r^ 
quires  velocities.  Assuming  that  axes  Xsj/s  rotates  with  the  angular  velocity  6s 
relative  to  the  inertial  axes,  the  velocity  vector  of  a  typical  point  in  s  can  be 
easily  derived.  Since  we  are  concerned  with  the  kinetic  energy,  let  us  derive  the 
velocity  squared.  Hence,  using  Eq.  (2.1),  we  can  obtain 

=  \wjw. 

=  i  jiJ  it,  -  Dj  cj  ii, 

(2-4) 
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where 


s9, 

—c9s 

■  0  -1  ■ 

D,  = 

c9. 

s9. 

,  1  = 

1  0 

in  which  c9  =  cos  6  and  s6  =  sin  9, 

For  chains  of  flexible  substructures  such  as  those  under  consideration  here, 
no  closed-form  solution  of  the  hybrid  differential  equations  is  possible,  so  that 
we  must  be  content  with  an  approximate  solution.  This  implies  discretization  in 
space  of  the  partial  differential  equations,  which  amounts  to  replacing  the  partial 
differential  equations  by  sets  of  ordinary  differential  equations.  The  discretiza¬ 
tion  can  be  carried  out  by  means  of  the  finite  element  method  of  the  classical 
Rayleigh-Ritz  method.  In  either  case,  we  express  the  elastic  displacement  vector 
in  the  following  form: 


Xs  0 

0 


(2.6) 


where  is  the  vector  of  admissible  functions  and  is  the  vector  of  gener¬ 

alized  coordinates.  We  will  assume  that  the  admissible  functions  introduced  in 
Eq.  (2.5)  automatically  satisfy  the  geometric  boundary  conditions.  In  addition, 
the  control  design  for  the  hybrid  equation  is  not  feasible,  so  that  the  use  of  the 
discretization  approach  can  be  justified. 

Inserting  Eq.  (2.5)  into  Eqs.  (2.4)  and  integrating  over  the  length  results  in 

T  =  ^m,RjR,-iRjDj{LS,+N^,q,)  +  RjCjN^,q, 

+ll,9^^+e,:$,q,  +  ^qjM,q,  (2.7) 

where 


m. 


r^t  rig 

=  /  rrisdxs  Ss  —  I  rfi^x^dXi  Is  =  nlsX^dxs  (2.8) 

Jo  Jo  Jo 


in  which  is  the  mass  density  per  unit  length.  In  addition, 

^s  =  f  rus^sdxs  $^  =  /  rfigXs^sdXs 

Jo  Jo 

M,  =  ^  =  [  J  ]  ^  =  J 


(2.9) 


(2.10) 


Moreover,  the  potential  energy  can  be  expressed  in  the  following  form; 


V,  = 


(2.11) 
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where  [  ]  denoted  symbplically  an  energy  integral  and  is  a  measure  of  the  strain 
energy.  Since  we  are  concerned  with  the  beam  bending  vibration,  the  elastic 
potential  takes  the  form 


d‘^Ui(x,t) 

dx] 


2 


dxt 


(2.12) 


Inserting  Eq.  (2.6)  into  Eq.  (2.12),  we  obtain  the  potential  energy  in  the  dis¬ 
cretized  form:  . 

(2.13) 

(2.14) 


where 


V,  =  -^qi  K,q, 

L  'dxl  dxl’^^ 


represents  the  substructure  stiffness  matrix. 

The  virtual  work  for  the  system  can  be  expressed  in  the  following  form: 


6W  =  T,6e, 


(2.15) 


where  T,  is  the  torque  applied  to  the  origin  of  the  body  axes  Xjj/j. 

Hence,  the  Lagrange’s  equations  of  motion  for  the  system  shown  can  be 


written  as 


dt  \dRj 

— =  r 
\dej  60  s  ' 

-I-  — =  0 

dt  \dqs)  dqs  dqs 


£ 

dt 


(2.16) 

(2.17) 

(2.18) 


Inserting  Eqs.  (2.7),  (2.13)  and  (2.15)  into  Eqs.  (2.16),  (2.17)  and  (2.18).  we 
can  obtain 


msRs  —  {LSs  + 

^Cj N'^sqs  -  (LSs  -h  N'^.qs)  -  20, =  0  (2.19) 

-  RJdKLS,  +  N^,q,)  +  Isl  +  =  T,  (2.20) 

CsR,  +  ^^0,  +  Msq,  +  K,q,  =  0  (2.21) 

The  above  equations  represent  the  nonlinear  ordinary  differential  equations  for 
the  single  substructure.  Imposing  the  kinematical  relations  to  the  above  equa¬ 
tions  to  link  two  adjacent  substructures  is  not  an  easy  task.  In  addition,  the 
control  design  for  the  nonlinear  system  is  not  easy.  Thus,  we  propose  the  per¬ 
turbation  method  in  the  next  section  to  ease  the  kinematical  synthesis  and 

numerical  calculations. 
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2.2  Perturbation  Method 

Designing  the  maneuvering  and  control  for  articulated  systems  of  substructures 
is  very  difficult,  especially  if  the  design  is  to  be  optimal  in  some  fashion.  The 
difficulty  can  be  traced  to  the  fact  that  the  system  is  nonlinear  and  of  high  or¬ 
der.  The  nonlinearity  can  be  attributed  to  the  rigid-body  motions  and  the  high 
order  to  the  elastic  motions.  The  perturbation  approach  is  based  on  the  simple 
observation  that  rigid-body  motions  tend  to  be  large  compared  to  the  elastic 
motions.  Consistent  with  this,  let  us  assume  that  the  translations  and  rotation 
of  the  body  can  be  divided  into  zero-order  terms  and  first-order  terms  in  mag¬ 
nitude.  Elastic  displacements  are  assumed  to  be  small  so  that  the  generalized 
displacements  associated  with  elastic  motion  can  be  regarded  as  a  first-order 
term.  Thus,  we  may  write 


Rs  =  R,o  +  Rsi  Os  =  ^.0  +  T,  =  r,o  +  (2.22) 

Inserting  Eq,  (2.22)  into  trigonometric  functions  yields  the  following  relations. 

c6s  =  -  sOsqOsi  sOs  =  sOsQ  -h  (2.23) 

which  lead  to 

Cs  =  Cao  ^  Dso6s\  Ds  =  Dso  +  CsqOsi  (2.24) 

Introducing  Eqs.  (2.22)  and  (2.24)  into  Eqs.  (2.19),  (2.20)  and  (2.21),  we  can 
obtain  the  zero-order  equations  of  motion 

f^sRso  ““  =  0  (2.25) 

-  S,L'^D,oR,o  +  I,e,o  =  Zo  (2.26) 

as  well  as  the  first-order  equations  of  motion 


"fT^sRsl  SgDgQLOsl  C^qN^sQs  ~  ^^soC^qLSsOsI 

-ie,oCjo  -  e%D,o)LS,0,i  -  (ZoDjo  +  0^,oCso)N^>q,  =  0  (2.27) 

—  S,L^D,oRsi  +  Is0si  +  —  SsL'^CjoRsO^si  “  —  Zi  (2.28) 


N'^CsoRsI  +  ^a^'al  +  —  *3^5  LJ D,oRsO0sl  + 

=  -^'^N^C,oR,o  (2.29) 

In  matrix  form,  we  have 

MsqXso  =  fsQ  +  dsQ  (2.30) 

+  CslXsl  +  fCslXsi  =  fgi  -f-  dsi  (2.31) 
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where 


XsO  = 


Rso 

OsO 


dso  = 


e%s,cToL 

0 


Rs\ 

■  0 

6si 

II 

< 

Tsi 

0 

r 

0 

1 

0 

(2.32) 

(2.33) 


(2.34) 


and 


fC,i  = 


M,o  = 


m,  -SsDjoL 

-S,L’^D,q  Is 


ms  -SsDjoL  &^s 

A4si  =  1  — Ssl^  Dso  I> 

C  N'^Cso  Ms 

’0  -2dfo  -2ioDjoN'^s 

C,i  =  0  0  0 

.  0  0  0 

0  -{esoCfo  -  e%Dso)LSs  -(OsoDJo  +  e%c^)N'^, 

0  -SsL'^CsoRso  -RsoDJqN^s 

0  -'^'^N^DsoRso  R> 


(2.35) 

(2.36) 


(2.37) 


(2.38) 


The  main  advantage  of  using  the  perturbation  method  is  in  that  the  zero- 
order  equations  can  be  solved  independently  of  the  first-order  equations.  The 
zero-order  equations  are  nonlinear  but  in  low  order.  Once  the  zero-order  equa¬ 
tions  are  solved,  the  solution  of  the  zero-order  equations  enters  into  the  first- 
order  equations  and  complete  the  numerical  calculations.  So  far,  we  are  con¬ 
cerned  with  the  single  substructure.  In  the  next  section,  we  will  see  how  to 
assemble  individual  equations  into  a  global  equations  of  motion  by  means  of  the 
kinematical  synthesis. 


2.3  Zero- Order  Kinematical  Synthesis 

The  derivation  of  the  equations  for  the  individual  substructures  was  carried  out 
in  Sec.  2.2  as  if  the  substructures  were  free  to  move  relative  to  one  another.  In 
reality,  however,  they  are  all  part  of  an  interacting  system  of  flexible  bodies, 
so  that  the  motions  of  the  various  substructures  are  coupled.  Indeed,  because 
the  various  substructures  are  hinged  to  one  another,  the  motion  of  a  given 
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hinge  is  accounted  for  several  times,  once  for  each  substructure  sharing  the 
hinge.  As  a  result,  the  motions  defined  by  the  independent  equations  of  motion 
derived  in  Sec.  2.2  contain  redundant  coordinates.  In  this  regard,  we  note  that 
hinges  constrain  the  translational  motions  but  leaves  the  rotational  motions  free. 
The  kinematical  synthesis  to  be  introduced  shortly  is  designed  to  eliminate  the 
surplus  coordinates.  We  propose  to  carry  out  the  kinematical  synthesis  for  the 
zero-order  and  first-order  equations  separately.  Then,  we  can  conclude  that 
the  zero-order  rigid-body  displacements  and  velocities  of  the  substructures  are 
related  as  follows: 


s-l 

1=1 

(2.39) 

5-1 

9so  =  ^10  +  A'+i,o 
»=1 

(2.40) 

The  first  and  second  derivatives  with  respect  to  time  can 

be  easily  obtained. 

5  —  1 

Rso  ”  ”  y  ^  D^QL£iOiQ 

*=i 

(2.41) 

5—1 

0sO  =  ^10  4-  /?>+!, 0 

X  =  1 

(2.42) 

5-1  5-1 

Rso  =  -Y1 

1=1  *=i 

(2.43) 

5—1 

OgQ  =  6io  -h  y^  A-}-i,o 

(2.44) 

i=l 


Let  us  confine  the  case  to  a  two-link  flexible  body  system  such  as  PACE. 
In  this  case,  the  origin  of  the  first  body  is  identical  to  the  origin  of  the  inertial 
axes  and  the  body  is  rotating  on  a  plane.  The  second  body  is  connected  to 
the  first  body  through  the  pivot  joint.  The  zero-order  system  amounts  to  the 
two-link  rigid  body  system.  Thus,  the  overall  degree-of-freedom  is  two.  The 
first  variable  can  be  chosen  as  an  absolute  angle  of  the  first  body  with  respect 
to  the  inertial  axes  and  the  second  variable  can  be  chosen  as  the  relative  angle 
between  the  body  axes  of  the  second  body  and  the  body  axes  of  the  first  body. 
Thus,  we  may  write  for  the  global  zero-order  displacement  as  follows: 

20  =  [^10  /?20]^  (2.45) 

The  acceleration  vector  of  each  body  can  be  generally  written  as 

x^o  ==  -£'a0^o  4"  dtfO  >  5=1,2  (2.46) 
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It  can  be  readily  seen  that 


Eiq  — 


O21  O21 

1  0 


for  the  first  body  and 


£*20  = 


-DiqIi£i  O21 
1  1 


(2.47) 


(2.48) 


for  the  second  body  where  Onm  represents  an  n  x  m  null  matrix.  Inserting 
Eqs.  (2.47)  and  (2.48)  into  Eq.  (2.30)  and  adding  the  results,  we  can  obtain  the 
global  zero-order  equation. 

=  do  (2.49) 

where 

2  2 

Mo  =  ^^eJqMsqEso  ,  do  =  —  Msodso)  (2.50) 

5=1  3=1 

It  can  be  observed  from  Eq.  (2.49)  that  the  numerical  simulation  of  the  global 
zero-order  equation  can  be  carried  out  independently  of  the  first-order  equation 
so  that  the  maneuvering  control  can  be  also  designed  independently.  Hence, 
the  maneuvering  motions  excite  beam  vibrations  which  in  turn  cause  perturbed 
motions  in  rotational  motions. 


2.4  First-Order  Kinematical  Synthesis 

T^he  kinematical  synthesis  for  the  first-order  equations  is  carried  out  in  the  same 
way  as  for  the  zero-order  equations.  The  derivation  of  the  constraint  relations 
is  tedious  and  is  omitted  here  for  brevity.  The  basic  recursive  relations  are 


Rsi  =  Z)  {CToNH^)qi  -  DjoLeAi]  (2.51) 

i=l 

a  — 1 

6si  =  ^11  +  ^2  (2.52) 

«=i 

The  first  and  second  derivatives  with  respect  to  time  are  as  follows: 


a-l 


i=l 


Rsi  =  X)  -  0ioDfoNM^i)Qi  -  DfoLiien  -  eioClLeAi 

3-1 

0sl  =  ^11  +  X 


(2.53) 

(2.54) 
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ijji 


Y,  -  noDfoN^i{ei)qi  -  eioDlNH^i)qi 


i=l 

-DfoL£i6ii  -  26ioCj^Uieii  - 


(2.55) 


5-1 

^5l  =  ^11  +  ^ 
t=l 


(2.56) 


If  we  confine  ourselves  to  the  two-link  flexible  body  system  as  we  did  in 
the  previous  section,  the  substructure  first-order  displacement  vector  can  be 
expressed  in  terms  of  the  global  first-order  displacement  vector. 


Xsi  =  EsiZi  ,  s  =  1,2 


(2.57) 


where 

Zl  —  [^11  q\  /?2192]^ 

It  can  be  readily  seen  that 


0 

K5 

02n 

O2I 

02n 

II 

1 

Oln 

0 

Oln 

1 - 

0 

3 

7nn 

Onl 

Onn 

for  the  first  body  and 


^21  = 


-DIM  Cfo/2^l(4)  O21  02n 

1  ^'l(^l)  0  0i„ 

0„l  Onn  Onl  Tfin 


(2.58) 


(2.59) 


(2.60) 


for  the  second  body  where  represents  an  n  x  n  identity  matrix.  Inserting 
Eqs.  (2.59)  and  (2.60)  into  Eq.  (2.31)  and  adding  the  results,  we  can  obtain  the 
global  first-order  equation. 


M\Zi  -h  CiZi  +  KiZ\  =  di  (2.61) 

where 

2  2 

Ml  =  (2.62) 

J  =  1  5  =  1 

2  2 
ICi  =  YeJi{M,iE,i  +  C,iE,i  +  }C,iEsi)  ,  di  =  (2.63) 

5=1  5=1 

It  can  be  observed  from  Eqs.  (2.61),  (2.62)  and  (2,63)  that  the  zeroorder  solu¬ 
tion  enters  into  the  coefficient  matrices  of  the  global  first-order  equation  so  that 
the  global  first-order  equation  is  in  fact  the  time-varying  linear  second-order 
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differential  equation.  Numerical  simulation  is  carried  out  in  a  way  that  first  we 
solve  the  global  zero-order  equation  and  then  insert  the  zero-order  solution  to 
the  coefficients  of  the  global  first-order  equation  and  solve  the  global  first-order 
equation.  It  should  be  noted  here  that  the  actual  rotational  motions  are  the 
sum  of  two  solutions  as  indicated  by  Eq.  (2.22). 
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Chapter  3 

Numerical  Results  and 
Discussions 


A  computer  program  based  on  the  algorithm  developed  here  was  written  using 
MATLAB.  We  consider  the  PACE  shown  in  Fig.  1.1  as  a  numerical  example. 
Each  body  represents  a  uniform  slender  beam  with  tip  masses  at  both  ends. 
The  maneuver  consists  of  bang-bang  controls  applied  to  each  DC  motors.  As  a 
result,  the  beams  undergo  bending  deformation  alone. 

Figures  3.1  through  3.4  show  time  histories  of  the  response  of  the  system. 
The  time  integration  of  equations  was  carried  out  by  means  of  the  wilson  6 
method.  Figure  3.1  shows  zero-order  response  of  the  system.  It  can  be  seen  from 
Fig.  3.1  that  the  rigid-body  rotations  are  made  through  the  torques  applied  to 
the  shoulder  and  the  elbow.  Figure  3.2  shows  first-order  rotational  responses  of 
the  system,  which  are  in  fact  the  per  tubed  motions  due  to  the  elastic  vibrations 
of  each  beam.  Thus,  the  actual  rotational  responses  of  the  system  are  the 
combination  of  Figs.  3.1  and  3.2  as  shown  in  Fig.  3.3. 

Figure  3.4  shows  time  responses  of  tip  displacements  of  each  beam.  As 
can  be  seen  from  Fig.  3.4,  the  fundamental  modes  are  mostly  excited  by  the 
maneuvering  motions.  However,  dynamic  characteristics  of  the  whole  system 
becomes  more  complicated  than  a  single-body  system  since  there  are  interactions 
between  two  bodies.  It  should  be  mentioned  here  that  an  efficient  numerical 
integration  scheme  should  be  developed  to  cope  with  the  nonlinearity  of  the 
system  equations. 

There  remains  a  question  regarding  what  kind  of  data  is  good  for  neural 
network  training.  Since  there  is  no  available  guideline  on  this,  we  generated  the 
data  set  by  applying  random  signals  to  the  shoulder  and  elbow  motors,  which 
results  in  random  motions.  It  should  be  mentioned  here  that  based  on  the  past 
PACE  experiments  the  motions  of  PACE  are  not  repeatable  which  may  cause 
serious  problems  in  neural  network  control  since  neural  networks  are  simply  the 
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\theta_{10}  and  \betaJ20} 


Figure  3.1:  Zero-Order  Displacements 


mapping  of  the  input  and  output.  When  there  exists  uncertainty,  there  is  no 
definite  answer  on  the  stability  of  neural  network  control.  We  will  have  ample 
opportunities  to  explore  the  neural  network  control  for  the  new  field  such  as 
control  of  flexible  multibody  systems. 

The  first  part  of  this  report  is  concerned  with  the  development  of  a  mathe¬ 
matical  formulation  capable  of  treating  the  problem  of  maneuvering  and  control 
of  flexible  multibody  systems.  The  formulation  is  based  on  equations  of  motion 
in  terms  of  coordinates  derived  for  each  substructure  independently  of  the  other 
substructures.  The  individual  substructures  are  made  to  act  as  a  single  structure 
by  means  of  a  consistent  kinematical  synthesis.  The  net  effect  of  this  synthe¬ 
sis  is  the  elimination  of  the  redundant  coordinates  resulting  from  the  original 
treatment  of  the  substructures  as  if  they  were  independent. 

In  most  problems  of  interest,  the  elastic  motions  are  small  compared  to 
the  maneuvering  motions.  In  recognition  of  this,  a  perturbation  approach  is 
developed  whereby  the  rigid-body  maneuvering  of  the  system  defines  the  zero- 
order  problem  and  the  elastic  motions  and  perturbations  from  the  rigid-body 
maneuvering  define  the  first-order  problem.,  where  the  term  “order”  refers  to 
magnitude  in  a  perturbation  sense.  The  kinematical  synthesis  mentioned  earlier 
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Figure  3.2:  First-Order  Displacements 


is  applied  to  both  the  zero-order  problem  and  the  first-order  problem.  The  zero- 
order  problem  for  the  rigid-body  maneuvering  is  nonlinear  and  of  relatively  low 
dimension  and  can  be  solved  independently  of  the  first-order  problem.  On  the 
other  hand,  the  zero-order  solution  induces  time- varying  coefficients  and  persis¬ 
tent  disturbances  in  the  first-order  problem.  The  system  of  equations  describing 
the  first-order  problem  is  linear,  time-varying,  and  of  high  dimension.  The  equa¬ 
tions  are  in  terms  of  components  about  body  axes,  which  makes  them  ideal  to 
control  design.  The  formulation  lends  itself  to  ready  computer  implementation. 

As  a  numerical  example,  PACE  test  article  is  maneuvered  according  to  a 
bang-bang  control  law  by  means  of  torques  applied  to  the  shoulder  and  the  el¬ 
bow.  The  numerical  example  shows  the  effectiveness  of  the  algorithm  developed 
in  this  report  as  well  as  the  applicability  of  the  equations  to  the  maneuvering 
and  control  of  flexible  multibody  systems. 
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and  tip_2 


of  Tip  Displacements 


Chapter  4 

Neural  Network  in  Control 


Contemporary  control  synthesis  algorithms  (e.g.,  robust,  adaptive,  and  optunal 
controls)  face  severe  limitations  for  some  of  the  more  challenging  realistic  systems  as 
summarized  in  (Ref.  22).  In  particular,  modem  space  structures  are  made  of  light¬ 
weight  composites  and  equipped  with  piezoelectric  or  piezoceramic  sensors  and 
actuators.  These  flexible  space  stmctures,  which  are  likely  to  be  highly  nonlinear  with 
time-varying  structural  parameters  and  poorly  modeled  dynamics  pose  serious 
difficulties  for  all  currently  advocated  methods.  The  difficulties  arise  in  the 
characteristics  of  control  system  designs  from  a  broad  spectmm  of  aerospace 
applications;  e.g.,  surveillance  satellites,  military  robots  or  space  vehicles.  Existing 
design  techniques  often  rely  on  the  assumption  of  a  good  dynamic  model  containing 
identified  system  parameters.  Furthermore,  these  iterative  and  computationally 
expensive  approaches  often  require  a  prior  fixed  design  constraints,  where  the  loading 
and  material  properties  need  to  be  specified  in  advances.  Consequently,  design 
procedures  to  achieve  the  desired  stability,  robusmess,  and  dynamic  response  for  large 
space  stmctures  with  unknown  parameters  are  incomplete. 

The  ultimate  autonomous  control  for  flexible  space  stmctures,  intended  to 
maintain  acceptable  performance  over  an  extended  operating  range  without  external 
intervention,  can  be  especially  difficult  to  achieve  due  to  factors  such  as  high 
dimensionality,  multiple  inputs  and  outputs,  complex  performance  criteria, 
operational  constraints,  imperfect  measurements,  as  well  as  the  inevitable  failures  of 
various  actuators,  sensors,  or  other  components.  It  is  impractical  to  expect  the 
existence  of  a  high  fidelity  model  prior  to  deployment  of  such  complex  stmctures. 
Therefore,  the  controller  needs  either  to  be  exceptionally  robust  or  adaptable  after 
deployment.  Also,  catastrophic  changes  to  the  stmctural  parameters,  due  to 
component  failures,  unpredictable  uncertainties,  and  environmental  influences, 
require  that  the  controller  be  reconfigurable  (i.e.,  it  refers  to  the  ability  of 
reconfiguring  the  controller  when  the  plant  dynamics  are  subjected  to  catastrophic 
changes.). 

Artificial  neural  networks,  based  on  simplified  models  of  the  human  brain,  have 
been  shown  with  promising  capability  for  data  manipulation.  They  have  proven  to 
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provide  an  efficient  means  of  learning  concepts  from  past  experience,  abstracting 
features  from  uncorrelated  data,  and  generalizing  solutions  from  unforeseen  inputs 
(Refs.  23-24).  Unlike  conventional  digital  computers  which  require  storing  program 
representing  the  processing  steps  need  to  be  followed  exactly,  artificial  neural 
networks  are  composed  of  highly  interconnected  processing  elements  which  are 
analogous  to  the  most  elementary  functions  of  the  biological  counterparts.  It  is  also 
believed  that  such  networks  are  extremely  robust  with  respect  to  malfunctions  of 
individual  devices  because  of  the  distributed  storage  and  information  flow.  They  are 
in  contrast  to  today’s  high  speed  Von-Neumann  computers,  which  heavily  rely  on  the 
perfect  functioning  of  each  device.  In  addition,  these  types  of  large-scale  dynamic 
systems  possess  inherent  potential  for  parallel  computation  because  of  assumed  highly 
interconnected  topology.  Indeed,  artificial  neural  networks  have  been  successfully 
applied  to  identification  and  control  of  many  dynamical  systems,  including 
autonomous  space  and  underwater  vehicles  (Refs.  25-26),  nuclear  power  plants  (Ref. 
27),  chemical  process  facilities  (Ref  28),  and  manufacturing  production  lines  (Ref 
29). 

The  ability  of  artificial  neural  networks  to  approximate  arbitrary  continuous 
functions  (e.g.,  nonlinear  plant  dynamics  and  complex  control  law)  provides  an 
efficient  mechanism  for  identification  and  control  of  large  space  structures.  Back- 
propagation  (BP)  networks  (Ref  30)  have  recently  been  demonstrated  to  have  the 
desired  functional  approximation  capability  with  an  arbitrary  degree  of  accuracy  (Ref 
31).  Although  back-propagation  learning  rule  proves  its  effectiveness  in  many 
instances,  it  usually  takes  excessive  time  to  train  the  neural  networks  and  the 
networks  may  get  trapped  into  local  minima. 

In  this  report,  we  utilize  an  efficient  paradigm  for  identification  and  control  of 
flexible  space  structures  by  incorporating  the  adaptive  time-delay  version  of  radial 
basis  function  (RBF)  networks  (Ref  32).  Radial  basis  function  network  with  one 
hidden  layer  has  been  shown  to  be  a  universal  approximator  (Ref  33).  These  model- 
free  neural  network  paradigms  (i.e.,  BP  and  RBF)  are  more  effective  in  solving 
control  problems  than  conventional  learning  control  approaches,  e.g.,  the  BOXES 
algorithm  which  partitions  the  control  law  in  the  form  of  a  look-up  table  (Ref.  34). 
Figure  4.1  shows  the  realizations  of  control  laws  via  the  BOXES  approach  (a 
staircase  approximation)  and  a  neural  network  approach  (a  smooth  approximation). 


Figure  4,1:  BOXES  and  Neural  Network  Approximations  of  Control  Law 
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Our  goal  is  to  approach  structural  autonomy  by  extending  the  control  system's 
operating  envelope,  which  has  traditionally  required  vast  memory  usage.  Artificial 
neural  networks,  on  the  other  hand,  overcome  common  memory  intensive  problems 
and  yet  provide  a  sufficiently  generalized  solution  space.  In  the  vibration  suppression 
problems,  we  utilize  the  adaptive  time-delay  radial  basis  function  network  as  a 
building  block  to  allow  the  control  system  to  function  as  an  indirect  closed-loop 
controller.  Prior  to  training  the  compensator,  a  neural  identifier  based  on  an  ARMA 
model  is  utilized  to  identify  the  open-loop  system.  The  horizon-of-one  predictive 
controller  regulates  the  dynamics  of  the  nonlinear  plant  to  follow  a  pnespecified 
reference  system  asymptotically  as  depicted  in  Fig.  4.2  (i.e.,  model  reference  adaptive 
control  architecture).  As  far  as  trajectory  slewing  problems  are  concerned,  the 
generalized  learning  controller  synthesized  by  the  adaptive  time-delay  radial  basis 
function  network  compensates  the  nonlinear  large  space  structure  in  a  closed-loop 
fashion  as  given  in  Fig.  4.3  (i.e.,  TDL:  tapped  delay  line). 


Figure  4.3:  Trajectory  Maneuvering  Learning  Control  Architecture 

The  function  of  the  neural  controller  is  to  map  the  states  of  the  system  into 
corresponding  control  actions  in  order  to  force  the  plant  dynamics  to  match  a  certain 
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output  behavior  which  is  specified  eitlier  by  the  reference  model  (Y,^)  or  command 
output  (Yjj).  However,  we  cannot  apply  the  optimization  procedure  (e.g.,  gradient 
descent  or  conjugate  gradient)  to  adjust  the  weights  of  neural  controller  because  the 
desired  outputs  for  neural  controller  are  not  available.  In  (Ref.  35),  a  specialized 
learning  algorithm  which  treats  the  plant  as  an  additional  unmodifiable  layer  of 
network  is  proposed.  The  output  error,  62,  is  back-propagated  through  the  plant  to 

derive  the  controller  output  error  .  and  can  be  approximated  by 

dYpj  Ypj(jl+Sujij)-Ypi(u) 

duj  5uj 

where  ^  is  the  ith  component  of  plant  output,  and  u-  is  the  )th  component  of  plant 

input.  However,  the  authors  fail  to  suggest  a  reliable  way  to  compute  ^2  •  In  (Ref.  36), 

the  inverse  Jacobian  of  the  plant  is  used  to  estimate  62  at  each  weight  update,  which 
results  in  a  complicated  and  computational  expensive  learning  procedure.  Moreover, 
since  the  plant  is  often  not  well-modeled  because  of  modeling  uncertainties,  the  exact 
partial  derivatives  cannot  be  determined.  In  (Ref.  37),  a  'dynamic  sign 
approximation’  is  utilized  to  determine  the  direction  of  the  error  surface,  assuming 
the  qualitative  knowledge  of  the  plant.  This  is  not  necessarily  the  case  in  space 
structure  applications  which  are  often  equipped  with  highly  correlated  parameters.  To 
achieve  the  true  gradient  descent  of  the  square  of  the  error,  we  use  ‘dynamic  back 
propagation’  (Ref.  38)  to  accurately  approximate  the  required  partial  derivatives.  A 
single-layer  adaptive  time-delay  RBF  network  is  first  trained  to  identify  the  open-loop 
system.  The  resulting  neural  identifier  then  serves  as  extended  unmodifiable  layers  to 
train  the  compensator  (i.e.,  another  single-layer  adaptive  time-delay  RBF  network).  If 
the  structural  dynamics  are  to  change  as  a  function  of  time,  the  back-up  neural 
identifier  would  require  the  learning  algorithm  to  periodically  update  the  network 
parameters  accordingly. 

In  Section  4.1,  back-propagation  network  is  first  outlined  providing  an  underlined 
issue  pertaining  to  the  learning  algorithm.  Advanced  from  these  derivations,  the 
interconnecting  topology  and  learning  algorithm  of  radial  basis  function  network  are 
introduced  in  Section  4.2.  This  is  followed  by  discussing  an  adaptive  time-delay 
radial  basis  function  network  in  Section  4.3  which  provides  an  effective  mechanism 
to  capture  most  of  the  spatiotemporal  interactions  among  the  structural  members. 
Simulation  study  of  PACE  test  article  is  performed  in  Section  4.4  to  demonstrate  the 
feasibihty  and  flexibility  of  the  present  architecture. 

4.1  Back-Propagation  Network 


A  given  back-propagation  network  (Ref  30)  can  be  completely  described  by  its 
interconnecting  topology,  neuronic  characteristics,  and  learning  rule.  The  individual 
processing  unit  performs  its  computations  based  only  on  local  information.  This 
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biologically  inspired  mechanism  provides  the  desired  collective  computational 
properties.  Netwoilc  with  sigmoidal  processing  units  has  recently  been  demonstrated 
to  have  the  desired  functional  approximation  capabilities  with  an  arbitrary  degree  of 
accuracy  (Ref.  31).  This  implies  that  back-propagation  network  is  dense  on  QD), 
where  C(D)  denotes  the  set  of  continuous  functions  in  the  closed  and  bounded  set  D. 
Based  on  this  mathematical  foundation,  back-propagation  networks  have  been 
successfully  applied  in  various  problems;  from  military  signaturing  to  medical 
diagnosis;  and  from  speech  recognition  and  synthesis  to  robot  and  autonomous 
vehicle  controls. 

The  ouQJut  of  the  ith  sigmoidal  neuron  in  the  fcth  layer  of  the  network  is  defined 
by 


\w^v^j-^+b^,  i  =  l,...,N* 

(4.2a) 

;=i 

vf=gf(H?).  i  =  l,....iV* 

(4.2b) 

where  gf :  (-«>,«>)  ->  (-1,1)  is  a  sigmoidal  function  (i.e.,  continuously  differentiable, 
monotonically  increasing,  and  g(0)=0),  A'*'"*  is  the  number  of  neurons  in  the  (A:-l)st 

layer,  W~  is  an  adjustable  parameter  representing  the  interconnection  weight  between 
the  output  of  the  yth  neuron  of  the  (A:-l)th  layer  to  the  input  of  the  <th  neuron  of  the 

Jfcth  layer  and  hf  is  the  external  bias  of  the  ith  neuron  of  the  Wi  layer.  Figure  4.4 
presents  a  functional  diagram  for  a  single  neuron. 


Figure  4.4:  A  Functional  Diagram  for  a  Single  Neuron 

The  interconnection  of  these  simple  processing  elements  into  a  multilayer  network 
given  in  Hg.  4.5  results  in  a  general  nonlinear  mapping.  Layer  1  is  referred  to  as  the 
input  layer  since  the  inputs  directly  effect  these  neurons.  The  last  layer  is  referred  to 
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as  the  output  layer,  since  the  output  of  the  network  is  the  output  of  these  neurons.  For 
a  single  layered  network,  each  processing  unit  is  an  output  unit;  thus  if  we  know  the 
desired  output  of  the  network  we  can  calculate  the  output  error  for  each  processing 
unit.  System  dynamics  for  single  layer  perceptron  network  are  given  by 

Oi  =  gi{Ui)  =  gii'ZWijIj+bi),  i  =  l,...,N^  (4.3) 

;=i 


where  Ij  is  the  input  to  the  network  and  is  the  number  of  input  neurons. 
Assuming  that  our  goal  is  to  minimize  the  sum  of  the  squares  of  the  output  errors 

1 

described  by  £  =— where  o,-  is  the  actual  ouq)ut  of  the  ith  processing 

2,=i 

unit  and  d,-  is  the  desired  output  of  the  ith  processing  unit,  then  gradient  descents  of  E 
with  respect  to  the  interconnection  weights,  w^,\<i<N^,\<j<N^,  and  external 

biases  h,,l<j  < are  achieved  by  adapting  each  parameter  according  to  equation 
(i.e.,  Tj  is  used  for  the  learning  rate.) 

Aw^=-Tl(d,.-o,)-^|^/,-.  (4.4a) 

Ah,.  =  -11(4-0,)^^  (4.4b) 

dUi 

output  layer 

hidden- 2  layer 

hidden- 1  layer 

input  layer 

O  linear  neuron  sigmoidal  neuron 

Figure  4.5:  A  Multi-layer  Back-propagation  Network 
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In  a  multilayer  network,  knowledge  of  a  desired  network  output  does  not  imply 
knowledge  of  the  output  error  of  each  processing  element  within  the  network.  The 
remaining  layers  of  the  network  are  referred  to  as  hidden  layers  since  the  outputs  of 
each  neurons  in  these  layers  are  hidden  from  the  outputs  of  the  network.  The  question 
is  how  to  adapt  the  weights  to  the  hidden  layers  when  the  desired  ouq)ut  of  the  hidden 
layer  is  unknown.  The  solution,  known  as  the  back-propagation  learning  algorithm  is 
outlined  below. 

As  done  for  a  single  layer  network  we  wish  to  adapt  the  interconnection  weights  to 
perform  a  gradient  descent  of  the  squared  error  at  the  output  of  the  network.  For 
processing  unit  in  the  output  layer  of  the  network,  the  weight  adaptation  equation  is 

given  by  Eq.  (4.4)  with  the  network  input  Ij  replaced  by  the  output  of  the  jth  neuron 
of  the  previous  layer 

=Ti^^vf  s  (4.5a) 

AbP  =  (4-5b) 

where  =  -(<j.  -Of)  .  The  introduction  of  the  5  parameter  will  simplify 

ouf 

the  presentation  of  the  algorithm  for  the  hidden  layers.  Each  can  be  thought  of  as 
the  output  error  backpropagated  through  the  network  to  the  input  of  the  /th  neuron  of 


the  ^  layer  (i.e.,  6^ 

cuy 


).  Consider  now  a  weight  which  does  not  connect  to  the 


output  layer  of  the  network,  (i.e.,  ,  where  l<k<p).  In  this  case  we  will  assume 

that  4^‘‘’4s  known.  From  Eq.  (4.5)  this  will  be  true  for  k  =  p-l.  Using  this 
information  as  a  starting  point,  the  algorithm  presented  below  can  be  used  to 
propagate  the  error,  in  an  efficient  manner,  backwards  through  the  network 

calculating  the  remaining  6*  *s.  The  gradient  of  the  square  output  error  with  respect  to 


is  described  by 


A  k_^E  _  p"’  dE 


du^  J 


dE  fC'  dE 


,5  duf  ^ 


A 

‘^dtr 


(4.6a) 


(4.6b) 


25 


where 


dE 


^4' 


■ik+l 


8Mf  dvf^  a«f  8Mf  ’ 

M.-1 

dbr 


(4.7a) 


(4.7b) 


and  is  the  number  of  neurons  in  the  (*+l)st  layer  of  the  network.  Thus,  the 
gradient  calculation  simplifies  to 


...»  8E  ,, 

iiNj^=n5,. 


<5*=| 

j=l 


Sg-(uf) 


8ttf 


(4.8a) 


(4.8b) 

(4.9) 


4.2  Gaussian  Function  Network 


A  radial  basis  function  network  (RBF)  (Ref.  32)  is  a  two  layer  neural  network  whose 
outputs  form  a  linear  combination  of  the  basis  functions  derived  from  the  hidden 
neurons.  The  basis  function  in  the  hidden  layer  produces  a  localized  response  to  input 
stimulus  as  do  locally-tuned  receptive  fields  in  our  nervous  systems  (Ref.  39). 
Gaussian  function  network,  a  realization  of  RBF  network  using  Gaussian  kernels,  is 
widely  used  in  pattern  classification  and  function  approximation.  The  output  of 
Gaussian  neuron  in  the  hidden  layer  is  defined  by 

M-=exp(— L- — ^),  i  =  l,...,iV‘  (4.10) 

2C7^^ 

where  u}  is  the  output  of  the  ith  neuron  in  the  hidden  layer,  x  =  ^ati,...,x^o]  is  the 

input  vector,  wj  =[h'|„...,h'^o]^  denotes  the  weighting  vector  for  the  ith  neuron  in 

the  hidden  layer  (i.e.,  the  center  of  the  jth  Gaussian  kernel),  a?  is  the  normalization 
parameter  of  the  ith  neuron  (i.e.,  the  width  of  the  yth  Gaussian  kernel),  and  and 
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are  the  numbers  of  neurons  in  the  input  and  hidden  layers,  respectively.  Equation 
(4.10)  produces  a  radially  symmetric  output  with  an  unique  maximum  at  the  center 
and  which  drops  off  rapidly  to  zero  at  large  radii.  That  is,  it  produces  a  significant 
nonzero  response  only  when  the  input  falls  within  a  small  localized  region  of  the 
input  space.  The  output  layer  equations  are  described  by 


y,  =  wfu\  i  =  l,...,iV^ 


(4.11) 


where  y,-  is  the  output  of  the  I'th  neuron  in  the  output  layer,  Wj  — 
denotes  the  weighting  vector  for  the  ith  neuron  in  the  output  layer,  u  =  , 


r 

r 


is  the  output  vector  from  the  hidden  layer,  and  is  the  number  of  linear  neurons  in 
the  output  layer.  The  output  layer  neurons  form  a  weighted  linear  combination  of  the 
ou^uts  from  the  hidden  layer. 

Gaussian  function  networks  have  the  property  of  undergoing  local  changes  during 
training,  unlike  back-propagation  networks  which  experience  global  weighting 
adjustments  due  to  the  characteristics  of  sigmoidal  function.  The  localized  influence 
of  each  Gaussian  neuron  allows  the  learning  system  to  refine  its  functional 
approximation  in  a  successive  and  efficient  manner.  Figure  4.6  provides  a  comparison 
of  back-propagation  network  (at  left)  and  RBF  network  (at  right)  in  a  pattern 
classification  problem.  Back-propagation  network  requires  at  least  4  hidden  and  2 
ou^ut  sigmoid  neurons,  while  Gaussian  function  network  needs  only  one  Gaussian 
hidden  neurons  and  2  ouq)ut  linear  neurons  to  correctly  classify  two  classes  of 
patterns  denoted  by  ‘x’  and  ‘o’,  respectively. 


Figure  4.6:  A  Comparison  of  Pattern  Classifications  Problem. 

The  network  effectively  positions  the  Gaussian  kernel  at  the  center  of  the  data,  and 
then  weights  and  thresholds  it  appropriately  to  produce  the  circular  decision 
boundary.  The  hybrid  learning  algorithm  (Ref.  32),  which  employs  the  K-means 
clustering  for  the  hidden  layer  and  the  LMS  for  the  output  layer  further  ensures  a 
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faster  convergence  and  often  leads  to  better  performance  and  generalization  compared 
to  back-propagation  learning  rule.  The  combination  of  locality  of  representation  and 
linearity  of  learning  offers  tremendous  computational  efficiency  in  real-time  adaptive 
control.  In  order  to  guarantee  the  generalization  capability,  a  cross  validation 
algorithm  is  used,  and  the  number  of  training  samples  is  ten  times  more  than  the 
theoretical  VCdim  bounds  (Ref  40).  K-means  algorithm  is  perhaps  the  most  widely 
known  clustering  algorithm  because  of  its  simplicity  and  its  ability  to  produce  good 

results.  The  normalization  parameters,  a?,  are  obtained  once  the  clustering  algorithm 
is  complete.  They  represent  a  measure  of  the  spread  of  the  data  associated  with  each 
cluster.  The  cluster  widths  are  then  determined  by  the  average  distance  between  the 
cluster  centers  and  the  training  samples, 

<7;^=—  X  i  =  (4.12) 

tfiixeQ”  “ 

where  0^  is  the  set  of  training  patterns  belonging  to  /th  cluster  and  is  the  number 
of  samples  in  0^-.  This  is  followed  by  applying  a  LMS  algorithm  to  adapt  the  weights 
in  output  layer.  The  training  set  consists  of  input/output  pairs  but  now  the  input 
patterns  are  processed  by  the  hidden  layer  before  being  presented  to  the  output  layer. 


4.3  Adaptive  Time-delay  RBF  Network 

Biological  studies  have  shown  that  variable  time-delays  do  occur  along  axons  due  to 
different  conduction  times  and  different  lengths  of  axonal  fibers.  In  addition, 
temporal  properties  such  as  temporal  decays  and  integration  occur  frequently  at 
synapses.  Inspired  by  this  observation,  the  time-delay  back-propagation  network  was 
proposed  in  (Ref.  41)  for  solving  the  phoneme  recognition  problem.  In  this 
architecture,  each  neuron  takes  into  account  not  only  the  current  information  from  all 
neurons  of  the  previous  layer,  but  also  a  certain  amount  of  past  information  from 
those  neurons  due  to  delay  on  the  interconnections.  However,  a  fixed  amount  of  time- 
delay  throughout  the  training  process  has  limited  the  usage,  possibly  due  to  the 
mismatch  of  the  temporal  location  in  the  input  patterns.  To  overcome  this  limitation, 
Lin  et,  aL  has  developed  an  adaptive  time-delay  back-propagation  network  (Ref.  42) 
to  better  accommodate  the  varying  temporal  sequences,  and  to  provide  more 
flexibility  for  optimization  tasks.  In  a  similar  spirit,  the  adaptive  time-delay  radial 
basis  function  network  is  proposed  in  this  study  to  take  full  advantages  of  temporal 
pattern  matching  and  leaming/recalling  speed. 

A  given  adaptive  time-delay  radial  basis  function  network  can  also  be  completely 
described  by  its  interconnecting  topology,  neuronic  characteristics,  temporal  delays, 
and  learning  rules.  The  individual  processing  unit  performs  its  computations  based 
only  on  local  information.  The  output  of  the  Gaussian  neuron  in  the  hidden  layer  is 
defined  by 
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Ix-wW 

«)  =exp(-”  ■#-).  j=U..,Ni  (4.13) 

2  0} 

where  u)  is  the  output  of  the  jth  neuron  in  the  hidden  layer  (denoted  by  the 

superscript  1),  x  is  the  input  vector,  w]  denotes  the  weighting  vector  for  the  ;th 

neuron  in  the  hidden  layer  (i.e.,  the  center  of  the  jth  Gaussian  kernel),  (Tf  is  the 
normalization  parameter  of  the  yth  neuron  (i.e.,  the  width  of  the  jth  Gaussian  kernel), 
and  is  the  number  of  neurons  in  the  hidden  layer.  Equation  (4.13)  produces  a 
radially  symmetric  output  with  a  unique  maximum  at  the  center  dropping  off  rapidly 
to  zero  for  large  radii.  Inspired  by  the  adaptive  time-delay  back-propagation  network, 
the  output  equation  of  ATRBF  networks  is  described  by 


y;(^n)  — SS 


where  wjij  denotes  the  weight  from  the  ith  neuron  in  the  hidden  layer  to  the  jth 
neuron  in  the  output  layer  with  the  independent  time-delay  is  the 


output  from  the  ith  neuron  in  the  hidden  layer  at  time  Lj^  denotes  the 

number  of  delay  connections  between  the  ith  neuron  in  the  hidden  layer  and  the  jth 
neuron  in  the  output  layer.  Shared  with  generic  radial  basis  function  networks, 
adaptive  time-delay  Gaussian  function  networks  have  the  property  of  undergoing 
local  changes  during  training,  unlike  adaptive  time-delay  back-propagation  networks 
which  experience  global  weighting  adjustments  due  to  the  characteristics  of 
sigmoidal  functions.  The  localized  influence  of  each  Gaussian  neuron  allows  the 
learning  system  to  refine  its  functional  approximation  in  a  successive  and  efficient 
manner.  The  hybrid  learning  algorithm  can  also  be  used  in  the  training  of  the 
adaptive  time-delay  radial  basis  function  network.  The  adaptation  of  the  ouq>ut 
weights  and  time  delays  are  derived  based  on  error  back-propagation  to  minimize  the 
cost  function. 


t  N  ^ 


where  djit„)  indicates  the  desired  value  of  the  jth  ouq)ut  neuron  at  time  t„.  The 
weights  and  time-delays  are  updated  step  by  step  proportional  to  the  opposite 
direction  of  the  error  gradient  respectively. 


29 


(4.16a) 

A4/=-n2-^T^.  (4.16b) 

where  Tjj  and  TI2  are  the  learning  rates.  The  mathematical  derivation  of  this  learning 
algorithm  is  straightforward.  The  learning  rule  can  be  summarized  as  follows. 


Aw?.,,  =  'r\i(dj{t„)-yj{t„))'Zu}(t„  -  t?;,,),  (4.17a) 

/=1 

At|,,  =-T!\2(djit„)-yj(t„))'Zwl  iu]  (t„  -t?.,).  (4.17b) 

/=1 


4.4  Simulation  Study 

The  5-second  training  data  for  system  identification  are  generated  directly  from  the 
modeling  analysis  in  response  to  uniform  random  inputs.  Due  to  the  time  limitation, 
at  the  current  stage  only  the  back-propagation  network  is  used  to  validate  the 
proposed  architecture  on  this  problem  domain.  Based  on  the  Pfs  experience  over  the 
past  two  years,  it  is  firmly  believed  that  with  adaptive  time-delay  radial  basis  function 
network,  the  network  will  converge  at  a  faster  rate,  provided  sufficient  number  of 
Gaussian  neurons  are  used  in  the  hidden  layer.  The  motivation  for  exploring  the 
back-propagation  network  in  this  study  is  to  serve  as  a  head  start  for  a  proper  choice 
of  the  number  of  Gaussian  neurons  in  the  ATRBF  network.  The  authors  will  continue 
to  pursue  in  this  issue. 

Two  sets  of  programs  in  MATLAB  and  gnu-C  are  developed  for  this  study.  They 
show  different  advantages  in  various  aspects.  The  flowchart  for  C-simulation  is  given 
in  Fig.  4.7.  The  dynamics  of  the  plant  are  assumed  unknown.  System  identification  is 
simulated  by  a  three  layer  back-propagation  network.  Two  hidden  layers  with  neurons 
15  and  10,  respectively  are  used  to  ensure  the  flexibility  to  approximate  arbitrary  non- 
convex  regions.  The  output  layer  with  linear  activation  function  is  used  to  scale-up  to 
arbitrary  range  for  the  state  variables  spans.  The  control  objective  is  to  achieve 
trajectory  slewing  as  well  as  vibration  suppression  along  the  motion. 
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Figure  4.7:  Simulation  flowchart  for  gnu-C  routines 


Our  strategy  is  to  relax  the  vibration  suppression  task  for  the  first  few  seconds.  The 
RMS  errors  are  accumulated  only  based  on  the  trajectory  following  components,  i.e., 
the  first  four  components.  Around  1.5  second,  we  begin  to  suppress  the  structural 
vibration  for  the  last  four  components  while  maintaining  the  progress  of  trajectory 
following.  However,  by  confining  the  desired  responses  to  strictly  zeros  will 
deteriorate  the  network  performance,  even  affecting  the  first  four  components.  The 
way  we  proposed  in  this  study  is  to  setup  an  implicit  exponential  delaying  envelop,  so 
the  network  will  smoothly  catch  up  the  requirements  as  we  expected.  Whenever  the 
ouqiut  falls  within  this  exponential  envelop,  we  will  have  the  zero  error,  meaning  no 
weight  adjustments  will  be  taken.  This  indeed  shortens  the  training  time  significantly. 
Trajectory  slewing/vibration  suppression  is  performed  by  another  back-propagation 
network  with  only  one  hidden  layer  (125  neurons).  The  output  layer  is  equipped  with 
hyperbolic  tangent  function  for  the  activation  function.  The  closed-loop  controller 
regulates  the  dynamics  of  the  PACE  structure  to  follow  the  desired  outputs  as  given 
below: 

6.(0=-^stoS,  (4.18a) 

(4.18b) 

and  tip  displacements  and  their  velocity  are  zeros,  of  course.  Figure  4.8  presents  the 
performance  of  the  neural  identifier  with  respect  to  all  output  variable  respectively,  in 
response  to  random  inputs  after  training  for  20000  trials.  It  took  roughly  a  day  to 
train  it  in  a  50-MHz  Pentium  computer.  Although  the  neural  identifier  learned  to 
match  the  open-loop  system  in  a  reasonable  time  frame,  the  compensator  took  more 
than  four  days  to  converge  to  reasonable  accuracy,  mean  square  error  0.000562. 
Figure  4.9  displays  the  closed-loop  performance  with  respect  to  all  output  variables, 
respectively,  in  response  to  an  impulse.  As  shown  in  the  figure,  the  neural  regulator 
has  learned  to  follow  the  specified  trajectory  and  then  damp  out  the  structural 
vibration  successfully. 
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Figure  4.8:  Open-loop  Response  of  PACE  test  article 


Chapter  5 
Conclusions 


This  paper  is  concerned  with  the  development  of  a  mathematical  formulation  capable 
of  treating  the  problem  of  maneuvering  and  control  of  flexible  multibody  systems. 
The  formulation  is  based  on  equations  of  motion  in  terms  of  quasi-coordinates 
derived  for  each  substructure  independently  of  the  other  substructures.  The  individual 
substructures  are  made  to  act  as  a  single  structure  by  means  of  a  consistent 
kinematical  synthesis.  The  net  effect  of  this  synthesis  is  the  elimination  of  the 
redundant  coordinates  resulting  from  the  original  treatment  of  the  substructures  as  if 
they  were  independent. 

In  most  problems  of  interest,  the  elastic  motions  are  small  compared  to  the 
maneuvering  motions.  In  recognition  of  this,  a  perturbation  approach  is  developed 
whereby  the  rigid-body  maneuvering  of  the  system  defines  the  zero-order  problem 
and  the  elastic  motions  and  perturbations  from  the  rigid-body  maneuvering  define  the 
first-order  problem,  where  the  term  "order"  refers  to  magnitude  in  a  perturbation 
sense.  The  kinematical  synthesis  mentioned  above  is  applied  to  both  th  zero-order 
problem  and  the  first-order  problem.  The  zero-order  problem  for  the  rigid-body 
maneuvering  is  nonlinear  and  of  relatively  low  dimension,  and  it  can  be  solved 
independently  of  the  first-order  problem.  The  control  can  be  open  loop  or  closed  loop. 
On  the  other  hand,  the  zero-order  solution  induces  time-varying  coefficients  and 
persistent  disturbances  in  the  first-order  problem.  The  system  of  state  equations 
describing  the  first-order  problem  is  linear,  time-varying  and  of  high  dimension.  The 
equations  are  in  terms  of  component  about  body  axes,  which  makes  them  ideally 
suited  for  control  design.  The  control  can  be  divided  into  two  parts,  one  closed-loop 
part  designed  to  control  transient  disturbances  and  another  part  designed  to  reject 
persistent  disturbances.  The  formulation  lends  itself  to  ready  computer 
implementation. 

As  a  numerical  example,  a  planar  tliree-beam  system  is  maneuvered  according  to 
a  bang-bang  control  law  by  means  of  a  torque  applied  to  the  center  beam.  Each  beam 
is  equipped  with  two  actuators  in  charge  of  suppressing  vibrations  and  perturbations 
in  rigid-body  motions;  they  exerted  forces  according  to  the  direct  velocity  feedback 
control  law.  The  numerical  example  shows  the  effectiveness  of  the  algorithm 
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developed  in  this  paper  as  well  as  the  applicability  of  the  equations  to  the 
maneuvering  and  control  of  flexible  multibody  systems. 

To  support  the  above  modeling  efforts,  an  effective  control  algorithm  needs  to  be 
design.  We  propose  the  use  of  a  neural  network-based  control  system  for  flexible 
multibody  systems.  The  proposed  neural  controller  is  designed  to  achieve  trajectory 
maneuvering  of  structural  member  as  well  as  vibration  suppression  for  precision 
pointing  of  space  structures.  To  this  end,  a  hybrid  connectionist  system  is  employed, 
which  consists  of  a  system  identification  neural  network  and  a  controller  neural 
network.  This  leads  to  a  robust  and  fault  tolerant  controller  readily  appUcable  to  space 
structures  of  interest  to  the  Air  Force. 

The  proposed  technology  can  directly  be  incorporated  into  various  structure 
control  designs  for  civilian  and  military  satellites.  In  a  similar  spirit,  the  proposed 
architecture  can  be  extended  to  dynamic  controls  of  space  engines,  underwater 
vehicles,  chemical  processes,  power  plants,  and  manufacturing  scheduling.  The 
proposed  technology  will  greatly  reduce  operational  and  developmental  costs,  and 
provide  high  performance,  fault  tolerant  control  systems  in  various  application 
domains. 
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